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ABSTRACT 



Context. While sub-micron- and micron-sized dust grains are generally well mixed with the gas phase in protoplanetary disks, larger 

grains will be partially decoupled and as a consequence have a different distribution from that of the gas. This has ramifications for 

predictions of the observability of protoplanetary disks, for which gas-only studies will provide an inaccurate picture. Specifically, 

criteria for gap opening in the presence of a planet have generally been studied for the gas phase, whereas the situation can be quite 

different in the dust layer once grains reach mm sizes, which is what will be observed by ALMA. 

Aims. We aim to investigate the formation and structure of a planetary gap in the dust layer of a protoplanetary disk with an embedded 

planet. 

Methods. We perform 3D, gas+dust SPH simulations of a protoplanetary disk with a planet on a fixed circular orbit at 40 AU to study 

the evolution of both the gas and dust distributions and densities in the disk. We run a series of simulations in which the planet mass 

and the dust grain size varies. 

Results. We show that the gap in the dust layer is more striking than in the gas phase and that it is deeper and wider for more massive 

planets as well as for larger grains. For a massive enough planet, we note that cm-sized grains remain inside the gap in corotation 

and that their population in the outer disk shows an asymmetric structure, a signature of disk-planet interactions even for a circular 

planetary orbit, which should be observable with ALMA. 

Key words, planetary systems: protoplanetary disks - hydrodynamics - methods: numerical 



1. Introduction 

While we und erstand the general scenario of planet formation 



via accretion ("Mizunol 
Wetherill 1990; Pollack et 



T98a 

5t alJ] 



IWetherilll \l9Sdt. iLissaueii [19871 
19961) . the devil is in the detail and 



the processes by which tiny sub-micron grains grow into plan- 
etesimals, the building blocks of planets, is not well understood. 
With over 400 extrasolar planets detected to date as listed in the 
extrasolar planet encyclopedical, we can start to do some mean- 
ingful statistics on t he types of pl anets that exist in our galaxy 
dAlibert et al.ll2005b iMarchilllOOTl) to help constrain theories of 
planet formation. 

A wealth of analytical as well as numerical studies of 
the formation of gaps by a planet in a gas disk have shown 
how the shape of the gap depends on properties of the planet 
(mass) as w ell as the gaseous disk (pressure scale height, vis- 
cosity). S ee 'Crida et al.' (200a) for the case of gap shapes and 
iPapaloiz ou et al. (2007) for a more general review on planet mi- 
gration and gap formation. 

The dust phase, however, has been shown to behave dif- 
ferently to the gas. Dust experiences a headwind from the 
pressure-supported sub-Keplerian gas and the induced drag 



force slows the dust and makes it settle to the midplane 
and migrate inwards. The magnitude of these effects depends 
strongly on the grain size and the disk density (Weidenschillin 
1977; Stemnski & Valaaeas 1996, 1997; Gar aud et al.1 l20d" 
iGaraud & Lia.2004:.Barriere-Fouchet et al...2005lK 
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In recent years, gap formation by a planet embedded in disks 
of gas and dust has been studied by several authors, for dif- 
ferent planet masses and different sizes of the solid particles. 
iPaardekooper & Mellemal (12004. 2006 a) determined the spatial 
distribution of 1 mm grains with 2D simulations in order to de- 
rive the smallest planet mass that would result in an observable 
gap with ALMA (Atacama Large Millimeter Array). In their 
work the dust was strongly coupled to the gas and responded 
indirectly to the planet gravity through the radial pressure gra- 
dients caused by a Neptu ne-mass planet in the gas disk that led 
to the formation of a gap. iMuto & Inutsukal (12009') investigated 
the effect of a low-mass planet on the dust distribution by inject- 
ing one dust grain at a tiin e and deri ved criteria for gap open- 
ing. ICieciel^g et aP (l2007h as well as iMarzari & Scholll (l2000h 
focused on alr eady formed planetesi mals while the gas phase 
is still present. ICieciel^g et al.l (l2007l) considered planetesimals 
down to 1 m in size. They studied the effect of spiral structures 
in the circumprimary gas disk triggered by a secondary compan- 
ion in a tight binary system in order to derive relative velocities 
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between planetesimals and determine whether planet formation 
is possible in such systems. 

The presence of planets, and the gap they create when they 
are massive enough, can also help constrain the global properties 
of the gas (temperature, density, viscosity) as well as those of the 
dust (grain size distribution, degree of settling) in the disks. This 
can be achieved by measuring t he width and brightnes s of the 
gap, ideally for each phase (Wolf & D'Angelo 2005: Crida et alJ 
l2006HFouchetetal.ll2007h . 

It has been shown that ALMA will be able to observe a plan- 
etary gap opened by a 1 Mj planet at 5.2 AU from a 1 Mq star at 
a distance of 140 pc under the naive assumption that gas and dust 
are well mixed (W olf et al...2002) . It will certainly be possible to 
observe other features related to the presence of the planet, such 
as warm dust in its vicinit y or even spiral waves for distances 
not exceeding 100 pc (Wolf & D'Ang elo 2005). Here we focus 
on the formation and features of the gap itself in the case of a 
massive protoplanet. 

Our previous simulations of dust evolution in a typi- 
cal C lassical T-Tauri Star (CTTS) disk (Barri ere-Fouchet et alJ 
I2005L hereafter BF05) showed that the thickness of the dust 
layer depends on grain size because different sized grains fall to 
the midplane at different rates. We distinguished three dynami- 
cal regimes for the dust: (1) almost uncoupled for large grains 
where the dust component follows slightly perturbed Keplerian 
orbits and keeps its 3D distribution (if initially 3D); (2) weakly 
coupled for intermediate-sized grains for which settling is very 
efficient; (3) strongly coupled for small grains where grains are 
forced to follow the gas motion. 

We have also investigated the formation of a gap by a planet 
immersed in a Minimum Mass Solar Nebula (MMSN) and 
showed that dust settling makes the gap much more striking in 
the dust layer than in the gas phase because of its reduced verti- 
cal extension for w eakly coupled grains (Maddison et al. 2007; 
iFouchet et al.ll2007l hereafter F07) . Indeed, t he criterion for ga p 
formation depends on the disk scale height (ICrida et al.ll2006h . 
Because of the size of particles in the weakly coupled regime for 
the particular case of the MMSN (1 m, see Sect. |2]l, that study 
had no direct application to observations. Images of dusty disks 
at infrared wavelengths do not probe such large particles, but in- 
stead trace the smaller, strongly coupled dust grains whose dis- 
tribution is similar to that of the gas. 

In this paper, we examine the gap formation in the dust layer 
of CTTS disks, which are spatially more extended and less dense 
than the theoretical MMSN case. Our ultimate goal is to use the 
results of our hydrodynamical simulations to produce synthetic 
images for ALMA. We consider in particular the effects of a 
massive planet in the outer cooler regions of the disk. Indeed, 
several planets at large distances from their star have been de- 
tected, such as th ose recently announced orbiting Fomalhaut 
(iKalas et al.ll2008l) . or HR 8799 (iMarois et al.ll2008l) . The extra- 
solar planet encyclopedia lists 10 planets with a semi-major axis 
larger than 20 AU, and 7 of them have minimum masses larger 
than 5 Mj. We study the dust distribution in the weakly coupled 
regime, which corresponds for these CTTS disks to a size range 
(100 yum to 1 cm) that can directly be probed by current and fu- 
ture (sub)millimetre instruments. 

The present paper is the first part of this work, presenting the 
hydrodynamical simulations. In Sect.|2l we discuss the gas-dust 
interaction in the presence of a planet. In Sect. [3] we describe 
the numerical method and simulation suite. Results are described 
in Sect. 31 while the analysis and explanations are presented in 
Sect.|5] We conclude in Sect.|6l In a forthcoming companion pa- 
per, we will use the simulations presented here to produce syn- 



thetic images for ALMA and present the most favorable observ- 
ing configurations to detect the gap and associated structures. 



2. Dust behaviour under the effect of gas drag 

Dust grains immersed in a gas disk experience a drag force 
which, in the Epstein regime ( valid fo r the nebula parameters 
and grain sizes we consider, see lBF05h . is given by 



F-D- — Av, 



where 



h 



PgCs 



(1) 



(2) 



and ma is the grain mass, fs the stopping time, Av = Vd - Vg the 
differential gas-dust velocity, pd the dust intrinsic density, s the 
grain size, pg the gas density and c^ th e gas sound speed. 

Early studies by Weidenschil Ungl (II977 ) showed that there 
exists a grain size for which radial migration is fastest, satisfying 
the condition fs^K - 1 where Qk is the Keplerian pulsation. This 
optimal grain size is thus given by 



^opt 



Pd^K 



(3) 



This depends on the nebula parameters, and in particular on 
the gas density profile. For a vertically isothermal disk in 
hydrostatic equilibrium, the gas density profile is po{r,z) = 
Pg(r, 0) e"'"^^''*''' , where the disk scale height is H(r) = 
Cs{r)/D.K{r). Integrating vertically yields the surface density 
2g(r) - y/lnpgir, 0) Hir). Consequently, the optimal size in the 
disk midplane is 



■^opt 



(r,0) 



yllnpi 



(4) 



This is typically ~1 m m- 1 cm for CTTS disks (iBFOSh and 
~1 m for MMSN disks dWeidenschillingll 19771: 'K)7). The three 
regimes mentioned in Sect. [1] correspond to grains much larger 
than iopt (almost decoupled from the gas), of comparable size 
to iopt (weakly coupled), and much smaller than .Sopt (strongly 
coupled). The intermediate regime, for which both radial migra- 
tion and vertical settling are efficient, is the most dynamically 
int eresting. 

iHaghighipour & BossI (l2003h have shown that solid particles 
drift towards pressure maxima. Indeed, at the inner edge of a 
pressure maximum, pressure increases with radius and the pres- 
sure gradient is positive, acting in the same direction as the stel- 
lar gravity. Gas has to flow super-Keplerian in order not to fall 
on the star It speeds up the dust that, with no internal pressure, 
would flow at Keplerian velocities. In order to conserve angular 
momentum, which increases with radius in a Keplerian disk, dust 
that is sped up and therefore gains angular momentum needs to 
migrate outwards. At the outer edge, this is all reversed and dust, 
slowed down by a sub-Keplerian gas, has to drift inwards. As a 
result, dust is concentrated in the pressure maximum. 

From those considerations, in a disk with an embed- 
ded planet we therefore expect dust grains to migrate to- 
wards the pressure maxima on either side of the gap and 
to accumulate the re. This behaviour has be e n observed by, 
e.g., Paardekoope r & MeUemal (|2004 l2006ah . iMaddison eFal] 
(.2007,) . and£Oi 
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3. Simulations 

3.1. Code description 

We have developed a 3D, two-ph ase (g as+dust) Smoothed 
Particles Hydrodynamic (SPH) code ('BFOSl). We use it to model 
a protoplanetary disk of mass Mdisk, orbiting a star of mass M*, 
and containing a planet of mass Mp on a fixed circular orbit of 
radius r^. The code units are chosen such that G - M*. - r^ - I. 
The disk is treated as vertically isothermal, implying that the 
cooling is very efficient. Any heat produced by viscous dissipa- 
tion or by stellar irradiation is radiated away much faster than the 
gas dynamical timescale. Disk self-gravity is not implemented: 
it is negligible for the low-mass disks we study. Gas and dust 
are treated as two inter-penetrating fluids coupled by gas drag in 
the Epstein regime. In each simulation we consider a population 
of uniform-sized grains which do not shatter or grow. For full 
details of the code, see BF05. 

We use the standard SPH viscosity (iMonaghanll 19891) with 
ttsPH =0.1 and ySspH = 0.5. The aspH term was introduced 
to remove sub sonic velocity oscillations that follow shocks 
dMonaghan & Gingold 1983) and the jSsph term damps high 
Mach number sh ocks and prevents partic le interpenetration 
(lMonaghanlll989l) . The q-sph was shown by iMonaghaiil (Il985h 
to produce shear and bulk viscosity. Our choice of artificial SPH 
viscosity terms leads to a Shakura-Sunayaev ass parameter for 
the viscosity of order 10"^. This is within the ra nge of values in- 
dicated by observations of protoplanetary disks dHartmann et al.l 
Il998t iKing et aD l2007 h. For a discussion on the val idity of 
these viscosity parameters in protoplanetary disks, see iFOTi In 
that previous work, we used cksph = 0.1 and /3sph = and 
showed that the resulting orss was only slightly smaller than with 
jSspH = 0.5. In this work, we consider a more massive planet, 
which will generate a stronger spiral density wave. As a result, 
we need to increase ySspH- 

SPH has clear advantages over grid-based codes, in particu- 
lar when investigating boundary-free problems such as 3D flared 
protoplanetary disks. However, within the simulation commu- 
nity, SPH is often considered to be a poor choice when mod- 
elling planet-disk interactions since the resolution decreases 
in the gap becaus e of th e reduced number of particles there. 
Ide Val-Borro et al] (12006 *) compared the results from grid-based 
and SPH codes when modelling the planet-disk interactions and 
found that SPH codes predict the correct shape of the gap, al- 
beit with less resolution in the low density regions and with 
weaker planetary wakes. It should b e noted, however, that the 
disk model of the grid-based codes of lde Val-Borro et aP (l2006l) 
used a rather low viscosity (which SPH cannot achieve). If one 
were instead to consider more turbulent disks, as expected in 
nature, that gap would be shallower and edges not as sharp in 
both SPH and grid-based codes. Authors modeling low-viscosity 
regions (e.g. dead zones) expect a gap with razor-sharp edges, 
while it is certainly not the case for the more viscous disks we 
consider. As mentioned above, we use o-ss - 10"^ as suggested 
by observations of disks. 

While the SPH community is substantially smaller than the 
grid-based community, people have been trying to improve the 
SPH technique by in cludin g high order algorithms in Lagrangian 
methods. Inutsuk^ (l2002h proposes a new formalism called 
GSPH, for Godunov Smooth Particle Hydrodynamics, which 
uses a Riemann solver to improve the treatment of shocks. One 
issue with this approach is that complicated equations of state 
are more difficult t o imple ment than when one relies on arti- 
ficial viscosity (see iMonagh an 1997). iMaron & HowesI (l2003h 
propose Gradient Particle Magnetohydrodynamics, but, as dis- 



Table 1. Simulation suite. 



Mp=0.1Mj 
Mp = 0.5 M, 
Mp = l Ml 

Mp = 5 Mj 



; = lOOyum s = I mm 



s = I cm 



cussed bv lPricd (|2004|) . it is not clear whether the increased com- 
plexity and computing cost is compensated by the gain in accu- 
racy. B0rve et al. (2009) propose Regularized Smooth Particle 
Hydrodynamics where the solution is mapped onto a regular 
grid. This approach reduces the numerical noise but increases 
the diffusivi t v. A n ovel technique has recently been presented 
by ISpringell (l2010l) . where a moving grid is set up by means 
of tessellation. This technique combines the advantages of both 
Eulerian and Lagrangian techniques and seems very promising. 
These developments are mostly proofs of concepts and need to 
be extensively tested. For this work we continue with the stan- 
dard SPH whose achievements and limitations have been clearly 
addressed in the literature. 



3.2. Simulation parameters 

We model a typical CTTS disk with M* = 1 M©, Mdisk = 
0.02 Mq and comprising 1 % dust by mass. The disk initially ex- 
tends from 0. 1 to 3 code units, which corresponds to 4-120 AU, 
with a planet at orbital radius rp = 40 AU = 1 code unit. The 
initial surface density profile, 2(r) - "Ep r^'', is taken to be flat 
(p = 0) to follow our previous work (lF07h . InitiaUy, Sopt is there- 
fore also uniform in radius, with a value of ~ 1 .5 cm in the mid- 
plane. The initial temperature profile is T(r) = Tor"*, with q = I. 
The disk aspect ratio is initially chosen to be H/r = 0.05 at rp. 
The intrinsic grain density, p^, is 1,000 kg m"^. 

As in F07, we embed the planet in a gas disk and evolve the 
system for 8 planetary orbits and then add the dust phase. The 
star is kept fixed at the origin. The evolution is then followed 
for a total of 104 planetary orbits. The disk is allowed to spread 
outwards and particles are only removed from the simulation if 
they go beyond 4 code units (160 AU). The inner radius is fixed 
at 0. 1 code units (4 AU) and particles crossing that limit are as- 
sumed to be accreted by the star. .Crida & Morbidelli(.2007.) have 
shown that the disk interior to the gap is dramatically depleted if 
the inner boundary condition is too close to the planet. Our test 
simulations show that an open inner boundary at 0. 1 rp avoids 
this depletion of the inner gas disk, as extensively discussed in 
F07. (see their Fig. 11). For full details of the initial setup, see 
lF07l All simulations contain 200,000 particles per phase. 

We run a series of simulations in which we vary the grain 
size s from 100//m to 1 cm, so that dust is in the weakly coupled 
regime, and the planet mass, Mp, from 0.1 to 5 Mj, in order to 
study how these parameters aff'ect the resulting dust distribution 
around the planetary gap. Here we investigate larger planet-to- 
star mass ratios than in F07, which produce a stronger spiral 
density wave. Table [Upresents the simulation suite. 



4. Results 

In this section, we describe the results of the simulations. 
Analysis and discussion of the results will be presented in the 
next section. 
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Fig. 1. Density maps in midplane (top) and 
meridian plane (bottom) cuts of the disks. The 
three leftmost columns show the dust density 
for s = 100 yum, 1 mm and 1 cm, from left 
to right, and the right hand column shows the 
gas density. The rows show simulations with 
M„ = 0.1, 0.5, 1 and 5 Mj, from top to bottom. 



In Fig. [T]we plot the dust and gas distributions in the xy and 
rz planes (where the center of frame is fixed on the star) after 
104 planet orbits, coloured by volume density for all simula- 
tions. Fig.|2]shows the azimuthally averaged radial surface den- 
sity profile and the midplane volume density of the dust for disks 
with a 5 Mj planet and grains sizes from 100 fim to 1 cm, and 



of the gas and dust for disks with s - I mm and planet masses 
from 0.1 to5 Mj. 

The gas is little affected by the dust phase (see Sect. l5.4l l and 
its distribution is similar for all grain sizes. Its density is shown 
in the rightmost column of Fig. [T] and the right panels of Fig. |2] 
for s - I mm and all planet masses. As expected (see Sect.lTJ, 
the more massive the planet, the deeper the gap. We find that a 
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Fig. 2. Azimuthally averaged surface density 
(top) and midplane volume density (bottom) 
profiles after 104 planetary orbits. Left: dust 
densities for Mp = 5 Mj and diiferent grain 
sizes. Center: dust densities for s = I mm and 
diiferent planet masses. Right: gas densities for 
s = I mm and diiferent planet masses. 



5 Mj planet carves a well-defined gap in the gas component of 
our CTTS disk, that planets with masses of 0.5 or 1 Mj create 
shallow gaps, and that a 0. 1 Mj planet only slightly perturbs the 
gas density profile. 

The dust phase, however, has a very different distribution de- 
pending on the grain size, as can be seen by comparing all sim- 
ulations for Mp = 5 Mj on the bottom rows of Fig. [T]and the left 
panels of Fig.|2] The disk extension is dramatically reduced as s 
increases from 100 fim to 1 cm: both its outer radius and vertical 
thickness decrease. A spiral density wave is visible on the face- 
on views for 100 //m and 1 mm grains (Fig. [1]), similar to that 
seen in the gas disk, while only a small section of it is visible 
in the much narrower disk of 1 cm grains. The gap is wider and 
deeper, with higher densities at its outer edge, for larger grain 
sizes, and in all cases more pronounced than in the gas phase. 
For 1 cm grains, because they are the most efficiently settled, 
the midplane volume density better shows their concentration at 
the gap outer edge, with a higher dust-to-gas ratio (~ 0.3) than 
the surface density does (~ 0.09, compared to the initial uni- 
form value of 0.01). The disk interior to the gap is virtually un- 
affected by the presence of the planet for 100 fim grains, but has 
almost disappeared for the 1 mm grains and is no longer present 
for 1 cm grains, for which a population of grains in corotation 
with the planet can be seen instead. In that latter case, the gap is 
slightly asymmetric and the outer disk appears eccentric (Fig.[T]). 

The same effect of grain size is seen for other planet masses 
in the three upper rows of Fig. [1] the dust disk's outer radius 
and vertical extension are smaller, whereas the planet's effect is 
larger, for 1 cm grains than for 1 mm grains. 

The effect of the planet mass on the dust phase for those two 
sizes can be observed in the two center columns of Fig. [T] and 
for 1 mm grains in the center panels of Fig. |2l The disk has the 
same radial extent but is slightly thinner for less massive planets. 
As is seen for the gas, the gap width and depth increase as Mp 
increases, and the gap is always more pronounced in the dust 
phase than in the gas. A very shallow gap is visible in the dust 
for a 0. 1 Mj planet, where only a slight perturbation is seen in the 
gas, and planets of 0.5 or 1 Mj akeady carve well-defined gaps 
in the dust disk. Contrary to the 5 Mj case, the disk interior to the 
gap is still present and particles in corotation are not observed. 



Figure [3] plots the time evolution of the azimuthally aver- 
aged surface density radial profiles for the three grain sizes for 
a 5 Mj planet. It shows that the disk's outer radius changes very 
little for 100 yum grains, decreases slowly for 1 mm grains, and 
very rapidly for 1 cm grains. The same behaviour can be seen 
for the density in the disk interior to the gap. For 1 cm grains, 
the displacement of the density peak towards the center indicates 
that the whole inner disk is accreted by the star The gap forma- 
tion shows the opposite trend as it is depleted rather rapidly for 
100 jum grains, more slowly for 1 mm grains, and the density 
peak at the planet's orbital radius indicating the 1 cm corotating 
grains does not evolve. 

Recently it has been shown that material can become trapped 
in the Lagrange points in disks with planets (e.g. in the dust- 
only+planet disks of Wolf etal. (2007) and in the gas+planet 
disks of de Val-Borro et alJ (1200 6)). Our simulations do not 
take the stellar wobble into account and thus cannot model 
Lagrangian points and produce any accumulation of matter 
there. Be that as it may, while gas is seen to be trapped at L4 
and L5 in the inviscid simulations of de Val-Borro et al. (2006), 
viscosity eventually removes gas from these points and indeed, 
their viscous runs do not show such features. We therefore do 
not expect to see accumulation of either gas or dust at L4 and L5 
in our viscous disks. 



5. Discussion 

5. 1 . Understanding the gap 

Most of the features described in Sect. |4] can be explained by 
comparing the grain size to the optimal size s^pt, whose radial 
profile in the midplane is shown in Fig.|4]at different times for 
Mp = 5 Mj. The closer their size is to ^opt, the more efficiently 
dust grains settle to the midplane and radially migrate. 

Initially at a uniform value of 1.5 cm, .Sopt decreases only 
slightly outside the gap and stays of the order of 1 cm out to 
~120 AU. Centimetre-sized grains therefore migrate and set- 
tle very rapidly, resulting in a compact and thin disk seen in 
Fig. [1] A similar behaviour is seen in the inner disk, which is 
very rapidly emptied as particles accrete onto the star, before Sopt 
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Fig. 3. Time evolution of the azimuthally av- 
eraged midplane radial volume density profiles 
for Mp = 5 Mj and different grain sizes. 
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Fig. 4. Evolution of the radial profile of the optimal grain size, 
iopt, in the midplane for Mp - 5 Mj. 



decreases significantly below 1 cm there (Figs.[3]and|4|. In the 
gap, however, Sopt drops very rapidly from its value of ~ 1 cm 
at dust injection (8 orbits) to below 1 mm in just 40 orbits, due 
to the rapid carving of the gap in the gas disk. Centimetre-sized 
grains are thus quickly decoupled from the gas and stay on their 
initial orbits in the gap, in corotation with the planet. Millimetre- 
sized grains settle and migrate less efficiently in the outer disk 
since their size is farther from .Sopt- In the inner disk, their dis- 
tribution evolves slowly at first, when .Sopt is still close to 1 cm, 
then their depletion accelerates as ^opt decreases and gets closer 
to their size, until only a tenuous annulus of dust is left at the 
end of the simulation. We anticipate it will disappear altogether 
at later times. In the gap, iopt drops rapidly to 1 mm and below, 
resulting in the fast migration of 1 mm-sized grains out of the 
gap. Finally, 100 yum grains are more strongly coupled to the gas 
and evolve slowly on either side of the gap, since their size stays 
well below ^opt at all times. The increase of the midplane den- 
sity in the outer disk reflects their moderate settling. In the gap, 
on the other hand, s^pt quickly approaches lOO/im, causing their 
efficient depletion. 

We note that the thickness of the dust layer is not vanish- 
ingly small despite the lack of explicit turbulence. Indeed, the 
final shape of the dust layer in these simulations, in both radial 
extent and thickness, depends on the grain size. This is due to 
a complex interplay between both settling and radial migration 
times and grain size. The disk can be thicker in places due to 
radial migration pushing more material inwards but a change in 
Sopt (due to density) further inwards can cause the dust to pile up. 
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Fig. 5. Number of SPH particles accreted by the star (left) and 
the planet (right) as a function of time, for Mp - 5 Mj. 



This dictates the shape of the dust disk during the disk evolution. 
As it is a highly non-linear process, it is difficult to predict the 
final shape of the dust layer 

Figure|5]shows the accretion of gas and of each grain popula- 
tion onto the star and the 5 Mj planet as a function of time. Both 
1 cm and 1 mm grains are accreted very efficiently onto the star, 
about twice as much as the gas, showing their strong inward mi- 
gration. The complete disappearance of the inner disk for 1 cm 
grains is confirmed by their accretion rate, which stalls after 
~ 80 orbits. The slightly lower accretion rate of 1 mm grains 
leaves a small population in the inner disk. The more strongly 
coupled 100 fj.m grains have an accretion rate comparable to that 
of the gas. Conversely, the accretion onto the planet is very weak 
for 1 cm grains, which are decoupled and stay in the corotation 
region on horseshoe orbits (see below), and for 1 mm grains, 
which quickly migrate out of the gap. The 100 fim grains take 
longer to leave the gap and consequently more are accreted by 
the planet, as is the gas which can still flow through the gap. 

The dust behaviour for other planet masses can also be un- 
derstood from similar ,Sopt profiles, which can easily be infeiTed 
from the gas surface density profiles in the top right panel of 
Fig.|2](vialog iopt (cm) ^ log S (kg m"^) - 1 .4 from Eq. Q). iopt 
does not change much in the outer disk, the shape and density 
profile of which are therefore the same for all planet masses. For 
Mp = 0.1 to 1 Mj, the gas density and consequently ,Sopt do not 
decrease much in the gap, preventing grains from decoupling and 
staying in the corotation region. Particles thus tend to evacuate 
the gap and move towards the pressure maxima at its edges. The 
inner disk is constantly replenished with grains migrating from 
the gap, slowing down its accretion onto the star and preventing 
its disappearance at the end of the simulations for 1 cm grains. 
Smaller grains have sizes smaller than .Sopt at all times interior to 
the gap, they therefore migrate very slowly and do not fall onto 
the star 

The 5 Mj planet triggers a strong spiral density wave in the 
gas (see Fig. [TJ which, incidentally, enhances the gas accretion 
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Fig. 6. Map of the dust-to-gas ratio, computed as the ratio of 
volume densities, in the disk midplane for Mp = 5 Mj and s - 
1 cm (logarithmic scale). 



onto the star fsee'Laug hlin et aljl997h . and is responsible for the 
lower gas density inside the gap compared to other planet masses 
(Fig. |2]). This spiral wave corresponds to a pressure maximum. 
Dust particles drift towards its location and the spiral pattern is 
also seen in the dust phase. The dust-to-gas ratio, shown in Fig.|6] 
for 1 cm grains in the disk midplane, is everywhere larger than its 
initial value of 0.01 and even reaches unity along the spiral wave. 
(It is still higher in the corotation region where there is little gas 
left.) This shows that the dust is even more concentrated by the 
ac tion of the d rag than the gas is in the spiral wave, as was found 
bv lRice et all (I2004D . 

The variation of the disk thic kness with planet mass c an be 
understood in light of the work bv lEdear & OuillenI (l2008h . who 
establish that there are strong vertical motions of the gas in the 
spiral arms, and they speculate that this stirs the dust vertically. 
These motions are naturally stronger for more massive planets, 
resulting in a thicker dust disk which we see in Fig.[I] Because 
the spiral arms are stronger closer to their launching point, we 
also expect vertical motions to be stronger close to the gap edge 
than in the outer disk and, as a result, the disk is thicker close 
to the gap for the larger planet masses. The less perturbed disks 
with lower mass planets have a regular flared shape that is thicker 
at large radii. 

The Lagrangian nature of SPH can be used to confirm the 
interpretation of the dust behaviour detailed above by follow- 
ing the trajectories of individual SPH particles. In the following 
discussion, we focus on the simulation with Mp = 5 Mj and 
s = 1 cm. We select a sample of dust particles initially in the 
upper layer of the disk at various radii and compute their in- 
stantaneous orbital elements (semi-major axis, a, eccentricity, 
e, and inclination, /) from their positions and velocities at each 
timestep. Figure|2]plots the time evolution of a, e, and / as well as 
the trajectories of these particles in the rz plane. Dots are plot- 
ted at regular time intervals (bottom panel of Fig. [T) to show 
how fast particles traverse different parts of their trajectories. 
The particles are initially on inclined elliptical orbits and settle 
to the midplane while they migrate radially. The efficient damp- 
ing of their inclination is clearly visible. Particles that started 
in the outer disk (PI, P2, P3, P4) migrate inwards and pile up 
at the outer gap edge (60-70 AU). Those particles initially be- 
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Fig. 7. Time evolution of orbital elements a, e, and / (top panels) 
and trajectories in the rz plane (bottom) of SPH dust particles 
initially at various radii in the upper layer of the disk, for the 
simulation with Mp = 5 Mj and s - I cm. 



tween that edge of the gap and the planet (at 40 AU) migrate 
outwards this time, towards the pressure maximum at the gap 
edge (P5, P6). Their eccentricities are damped during their mi- 
gration, but when they reach the gap edge (which the outermost 
PI does not achieve by the end of the simulation), e increases 
again and oscillates about an average value between 0.04 and 
0.08. The eccentricities of P4 and P6 behave very similarly to 
that of P5 but are not plotted to avoid overcrowding the figure. 
The trajectory of P7, also shown in the xy plane in Fig.[8l is rep- 
resentative of dust grains that started in the vicinity of the planet: 
they are in corotation and follow horseshoe orbits. Finally, par- 
ticles initially interior to the corotation region migrate inwards 
towards the inner gap edge and end up being accreted by the star 
(P8, P9). Their eccentricities, also not plotted, oscillate between 
and 0.05. 

The positions of SPH dust particles at dust injection and at 
the end of the simulation are plotted in Fig. |9] Selecting those 
that are in the corotation region (in green) at the end of the simu- 
lation and looking for their initial positions show that they were 
akeady there. This confirms that they decouple very quickly 
from the gas and do not leave the gap, and that particles initially 
elsewhere do not migrate into the corotation region. 

Some particles initially in the corotation region can pass 
close to the planet, which then scatters them away. The trajec- 
tories of some of these particles up to their ejection are plotted 
in Fig. [8] (the coarseness of the curves is due to the limited fre- 
quency of the code data output), while Fig. [10] shows the time 
evolution of their orbital elements. PIO is initially on an inclined 
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Fig. 8. Trajectories in the xy plane, in the planet's reference 
frame, of SPH dust particles initially in the corotation region, for 
the simulation with Mp = 5 Mj and 5=1 cm. The Saturn-like 
symbol shows the planet's location. 




Fig. 9. Positions in the xy (top) and rz (bottom) planes of the 
SPH dust particles at dust injection (left) and at the end of the 
simulation (right) for Mp = 5 Mj and s -\ cm. Particles ending 
up in corotation are highlighted in green. 



orbit and settles to the midplane and keeps a horseshoe orbit for 
some time, then passes close to the planet and is scattered in- 
wards. By that time, there is not enough gas left in the inner disk 
to make it migrate further (iopt has decreased below 1 cm and the 
particle is decoupled from the gas) and it stays in a marginally 
perturbed orbit around 30 AU. Pll passes close to the planet 
early in the evolution, before completing a full horseshoe or- 
bit, and is scattered inwards. It keeps migrating to the inner gap 
edge, slightly outside 10 AU, where it stays for some time before 
eventually being accreted by the star together with the rest of the 
inner disk particles. P12 is scattered outwards with a high eccen- 
tricity after a long period on its horseshoe orbit. It then migrates 
inwards, its semi-major axis and eccentricity being damped by 
the gas. 
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Fig. 10. Time evolution of orbital elements of SPH dust particles 
initially in the corotation region, for the simulation with Mp = 
5 Mj and s -\ cm. 
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Fig. 11. Disk eccentricity as a function of radius for the gas and 
dust phases for Mp = 5 Mj and s - \ cm. 



5.2. Understanding the disk asymmetry 

Figures [1] and |9] show that the gap and outer disk become asym- 
metric in the dust phase for 1 cm-sized grains with a 5 Mj planet. 
As can clearly be seen in the gas phase of Fig. [1] the more mas- 
sive the planet the stronger the resulting spiral wave. This wave 
is not a "physical" wave in that it comprises specific particles 
throughout the simulation, but is instead a density wave which 
moves through the disk with a pattern speed that matches that 
of the planet's orbit. Thus it is not coherent eccentric orbits of 
the individual particles, but the density wave itself which creates 
the asymmetric distribution in the dust phase. The 1 cm particles 
respond most strongly to the pressure maximum in the density 
wave and therefore we see a stronger disk asymmetry. This is 
enhanced by the fact that the disk of 1 cm-sized grains is more 
co mpact, making th e stru cture more obvious. 

iKIev & DirksenI (l2006h studied gas-only disks with an em- 
bedded planet on a circular orbit and observed a transition from a 
circular to an eccentric disk for planet masses larger than 5.25 Mj 
in an o-ss = 0.01 disk. They found that the disk is turned eccen- 
tric by the excitation of an m - 2 spiral wave at the outer 1:3 
Lindblad resonance. The asymmetry in our dust disk cannot have 
the same origin because our planet masses are smaller than that 
necessary for an eccentric disk, and because our simulations are 
run for 100 orbits, which is much shorter than the viscous time 
scale necessary for eccentricity growth. We computed the disk 
eccentricity, defined as the average of the inst antaneous eccen- 
tricity of all particles at a given radius and which lKlev & DirksenI 
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Fig. 12. Torques due to the gas and dust phases on the planet in 
the simulation with Mp = 5 Mj and s = 1 cm. Left: normalized 
torques as a function of the radius in the disk at the end of the 
simulation. The planet's orbit is marked by the dotted line. Right: 
normalized total torque on the planet as a function of time. 



(2006") used as a diagnosis tool, and plot it for gas and dust in 
Fig. [TT] For both phases, the disk eccentricity reaches a maxi- 
mum of only about 0.08, similarly to their non-eccentric disks 
(see their Fig. 2). Finally, we do not see any asymmetry in our 
gas disk. The one we observe in the dust disk is therefore spe- 
cific to the dust dynamics in the presence of a spiral wave in the 
gas. 

As expected from the analytical study of iMarzari & Scholll 
(1200 0*). we do not observe periastron alignment of our dust par- 
ticles on neighbourin g orbits because our pl anet is not eccentric. 
However, contrary to lCieciel^g et al.l (l2007h we do not see a cor- 
relation between periastron longitude and eccentricity. This is 
likely due to the fact that our disk perturbation is much weaker 
due to the lower mass of our secondary compared to theirs. 

The disk asymmetry has implications for observations. 
iTakeuchi & Artvmowiczl(l2001h suggested that consistent detec- 
tions of planets through structures in the dust layer would be 
more convincing if asymmetric. They showed that it was possi- 
ble to build axisymmetric dusty rings and gaps without the pres- 
ence of a planet. Here we show that even a planet on a circular 
orbit can lead to an asymmetric dust structure, and for a grain 
size that will be observable by ALMA. 

5.3. Understanding torques 



In IF07I we suggested that the settled dust may have a non- 
negligible effect on planet migration. While the dust content is 
only 1 % in mass of that of the gas, once it has settled to the mid- 
plane where the planet orbits, it could exert a rather strong torque 
on it. We therefore computed torques exerted by the dust and 
compared it to those exerted by the gas in Fig. [T2]for our simu- 
lation with Mp = 5 Mj and s = I cm. In the left panel, we plot 
the torques as a function of radius at the end of the simulation. 
The torque exerted by the gas shows a positive component inside 
the orbit of the planet and a negative one outside, a well-known 
feature of planet gaps. The oscillations at larger radii are due 
to the spiral wave. The inner component is substantially weaker 
than the outer one because the inner disk is partly depleted. The 
torque exerted by the dust phase has the same global shape ex- 
cept that the inner component vanishes because the inner dust 
disk is empty. It is an order of magnitude smaller than that due 
to the gas because the densities also differ by an order of magni- 
tude. For 1 mm and 100 yum grains, the dust-to-gas ratio is even 
smaller, resulting in an even weaker dust torque (not shown). We 
therefore do not expect dust structures to have a direct dynamical 
influence on planet migration. 




Fig. 13. Azimuthally averaged midplane density profiles after 
104 planetary orbits for Mp = 1 Mj and 5=1 mm (top) or 
Mp - 5 Mj and i = 1 cm (bottom), for simulations with or 
without backreaction (BR) of dust on gas. Left: dust; right: gas. 



In the right panel of Fig. [12] we plot the total torque as a 
function of time. Here we can see that the absolute value of the 
torque exerted by the gas steadily decreases while that of the 
dust slowly increases. This is expected given that the gap in the 
gas phase is being carved until it reaches a steady state after 
several hundreds of orbits. Thus, the gas in the Lindblad res- 
onances (which are inside the Hill radius) is depleted and the 
Lindblad torque decreases. The dust, on the other hand, is pil- 
ing up at the outer gap edge, which results in the increase of 
the outer Lindblad torque. Eventually, as the system reaches a 
steady state, the increase in the Lindblad torque exerted by the 
dust phase will stop as well. Perhaps, when the stationary case 
is reached, will dust have a direct dynamical influence on the 
planet. 

5.4. Effect of backreaction of dust on gas 

The simulations presented so far are all run with no backreaction 
of the dust phase on the gas. To see whether this simplification 
is valid, we plot in Fig.[T3]the gas and dust volume densities for 
simulations with and without backreaction of dust on the gas, for 
two configurations: (Mp = 1 Mj, i = 1 mm) and (Mp = 5 Mj, s = 
1 cm). 

In the (1 Mj,l mm) case, the eff'ect of dust on the gas is 
mostly visible in the inner disk, where it drags the gas along with 
it towards the central star, decreasing the gas density. Indeed, af- 
ter 104 orbits, the number of gas SPH particles accreted by the 
star is about 10% larger with than without backreaction. In the 
dust, when backreaction is included, we see that the depth of the 
gap and the height of the outer edge are smaller and the outer 
disk radius is larger. This is because the backreaction means that 
the dust tends to drag the gas along with it, which effectively re- 
duces the gas drag on the dust. Therefore, its effects on the dust 
distribution (radial migration, gap creation, pileup at the outer 
edge) which we have discussed extensively are weakened. This 
means that the dust depletion in the gap is not as strong, the dust 
pileup at the gap edges is weakened and the radial migration in 
the outer disk is retarded. 

The effects for the (5 Mj,l cm) case are the same, and one 
may expect them to be more important for a larger planet mass 
causing stronger pressure gradients and for the grain size caus- 
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ing the largest dust-to-gas ratio. However, they are harder to see, 
for several reasons. The effect on the gas in the inner disk is 
smaller because all the dust there has been accreted onto the 
star well before the end of the simulation, and the time during 
which the backreaction of dust on gas can act is much smaller 
than in the (1 Mj,l mm) case. The effect of backreaction on the 
dust in the inner disk cannot be evaluated since the dust has dis- 
appeared there. In the gap, dust grains decouple from the gas 
very early on and stay on horseshoe orbits, their distribution is 
no longer affected by gas drag, with or without backreaction, and 
both curves are very similar With backreaction, the height of the 
outer edge in the dust is slightly reduced, as for the (1 Mj,l mm) 
case. Finally, the dust disk outer radius is also slightly larger with 
backreaction, but this is harder to see due to the steeper profile. 
The overall effects of backreaction in the simulations pre- 
sented in this paper are quite small, however, which validates 
our assumption of no backreaction for the simulations presented 
in the rest of this study. We do not investigate in this work the 
potential development of the streaming instability, which, in the 
presence of backreaction of dust on gas, can lead to density en- 
hancements and clumping in the dust layer (for more details, see 
fYoudin & Goodmanl2005MYoudin & Johansenll2007l) . 

5.5. Putting our work in context 

Over the past few years there has been a wealth of literature on 
the effects an embedded planet has on the gas and dust grains 
of various sizes in protoplanetary disks. lEdgar & OuillenI (|2008|) 
studied the vertical structure of spiral arms raised by an em- 
bedded planet, suggesting that the strong gas vertical motions 
should make the dust layer slightly thicker close to the gap (or at 
least over the radial extent of the spiral wave). Our simulations 
of gas and dust disks confirm their expectations as discussed in 
Sect. 15.11 which can have implications for the visibility of the 
gap- 



iPaardekooper & Mellemal (|2004 l2006ah studied the effects 
of low mass planets (Mp - 0.01 to 0.5 Mj) on strongly cou- 
pled grains (1 mm in size but with a denser and more compact 
nebula than in our study) in vertically integrated 2D disks. They 
were therefore unable to follow the vertical settling of dust and 
to investigate the variation of the disk thickness with grain size 
and the interplay between vertical and radial migration as we do. 
Because their grains were strongly coupled to the gas, the struc- 
tures in their disks where mostly due to the planet perturbation 
on the gas. Decoupling between gas and dust happened at shock 
fronts in the gas where it was decelerated while dust, with no 
internal pressure, was not. As a result, the flux of dust grains in- 
side the spiral arms was larger than that of the gas and, because 
spiral arms deflect particles away from the planet, the gap in the 
dust layer became deeper than that in the gas phase. In our case, 
planets are massive enough to directly deflect particles, which 
subsequently stay out of the gap (see Sect. 15.1b . Indeed, their 
eccentricity acquired during the close encounter with the planet 
is rapidly damped by the gas and they slowly migrate towards 
the outer gap edge, if deflected outwards (like particle PI 2 in 
Fig. [Jo]). 

Additionally, while IPaardekooper & MeUemal ( |2004 l2006al) 
observe resonant trapping (throug h an i ndirect mechanism), we 
do not (this has been discussed in lF07h . This difference can be 
traced back to the fact that our dust is marginally decoupled 
while theirs is strongly coupled despite the fact that we both 
study 1 mm grains (again our nebula parameters differ). As a re- 
sult, they expect multi-ringed structures which could be observa- 
tionally difficult to disentangle from a multi-planet system. We, 



on the other hand, see structures that can only be related to a sin- 
gle planet. We expect that smaller grains would behave similarly 
to what was found by Paardekooper & Mellema (2004, 2006^ 
but, in the submillimetre wavelength range, those smaller grains 
do not contribute much. We thus conclude that, for the nebula 
parameters we chose, there will be no ambiguity between single 
or multiple planet systems. 

Finally, as well as potentially exerting a torque on the planet, 
changes in the dust density will affect the dust opacity, which 
strongly influences the temperature structure around the planet. 
We study grains in a size range where their opacity is not dom- 
inant, but due to collisional shattering, they most likely give 
rise to a population of smaller grains with large opacities who 
will tend, at least to first approximation, to follow the spa- 
tial distribution of their larger parent grains. Recent studies us- 
ing radiative transfer in the flux-limited diffusion approximation 
(IPaardekooper & MeUemall2006bi l2008l) show that the tempera- 
ture structure in the disk and specifically around the planet has a 
strong e ffect on the migration rate, bu t also on the migration di- 
rection. iHasegawa & Pudnta ( 1201 Obi) have shown that dust set- 
tling and the presence of a dead zone have an effect on the tem- 
perature distribution. They note the formation of a dust wall at 
the edge of the dead zone that is directly illuminated by the star 
and see a local positive temperature gradient in front of the dust 
wall. They demonstrate that it has an effect on the migration 
of sma ll planets in a forthcoming paper (Hasegawa & Pudritg 
1201 Oah . It is however beyond the scope of our paper to study the 
variations in opacity and subsequent effect on planet dynamics. 



6. Conclusions 

We have run 3D, two-phase (gas and dust) SPH simulations of 
a typical CTTS disk of mass 0.02 Mq with an embedded gi- 
ant planet on a circular orbit at 40 AU. We vary the grain size 
(100 yum, 1 mm, 1 cm) and the planet mass (0.1, 0.5, 1, 5 Mj) 
and study the formation of the planetary gap. We confirm that 
gap opening is stronger in the settled dust layer than in the flared 
gas disk. Gaps are deeper and wider for (1) larger, more effi- 
ciently settled grains and (2) more massive planets. Larger planet 
masses are required to open a gap in the gas phase than in the 
dust: while a 0.5 Mj planet only slightly affects the gas phase, it 
carves a deep gap in the dust. For the most massive 5 Mj planet, 
1 cm grains remain trapped in corotation with the planet while 
their distribution in the outer disk shows an asymmetric struc- 
ture, even though the planet's orbit is circular. We find that this 
is not caused by the periastron alignment of a coherent set of or- 
bits but rather by the pile-up of dust in the pressure maximum of 
the gas phase caused by the spiral density wave triggered by the 
planet. This global asymmetry does not appear for less massive 
planets because the spiral perturbation is substantially weaker. 

The variety of structures that we obtain in the dust phase for 
various grain sizes and planet masses has implications on the ap- 
pearance of protoplanetary disks at (sub)millimetre wavelengths 
and show how important it i s to go beyond the gas-onl y disk de- 
scription proposed by, e.g., IWolf & D' Angelor(l2005l) . The ob- 
servability of these disks with ALMA is the subject of a forth- 
coming companion paper. 
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